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Abstract 

Multidimensional optimization problems where the objective function 
and the constraints are multiextremal non-differentiable Lipschitz functions 
(with unknown Lipschitz constants) and the feasible region is a finite col- 
lection of robust nonconvex subregions are considered. Both the objective 
function and the constraints may be partially defined. To solve such prob- 
lems an algorithm is proposed, that uses Peano space-filling curves and the 
index scheme to reduce the original problem to a Holder one- dimensional 
one. Local tuning on the behaviour of the objective function and constraints 
is used during the work of the global optimization procedure in order to accel- 
erate the search. The method neither uses penalty coefficients nor additional 
variables. Convergence conditions are established. Numerical experiments 
confirm the good performance of the technique. 
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1 Introduction 



In the last decades there has been a growing interest in approaching Global Opti- 
mization problems by different numerical techniques (see, for example, Hl- lfTTTl . 
II2I-II2OI, 1251 1211, Ea-Eg and references given therein). Such an interest is 
motivated by a large number of real-life applications where such problems arise 
(see, for example, [1H31, 10 El [8] [JOl [B] III] |22l |22 |25l |28l [30l These 
problems often lead to deal with multiextremal non-differentiable objective func- 
tion and constraints. In such a context the Lipschitz condition becomes the unique 
information about the problem. 

It has been proved by Stephens and Baritompa [37] that if the only informa- 
tion about the objective function (p{y) is that it belongs to the class of Lipschitz 
functions and the Lipschitz constant is unknown, there does not exist any deter- 
ministic or stochastic algorithm that, after a finite number of evaluations of ip{y), 
is able to provide an underestimate of the global minimum ip{y*). Of course, this 
result is very discouraging because usually in practice it is difficult to know the 
constant. Nevertheless, the necessity to solve practical problems remains. That is 
why in such problems instead of the statement (SI) 'Find an algorithm able to stop 
in a given time and provide an e-approximation of the global optimum f{y*)' the 
statement (S2) 'Find an algorithm able to stop in a given time and return the low- 
est value of ip{y) obtained by the algorithm' is used. Under this statement, either 
the computed solution (possibly improved by a local analysis) is accepted by final 
users (engineers, physicists, chemists, etc.) or the global search is repeated with 
changed parameters of the algorithm. 

Theoretical analysis of algorithms (depending on parameters) for solving prob- 
lems (S2) is similar to analysis of penalty methods. It is proved for global opti- 
mization algorithms that for a fixed problem there exists a parameter P* such that 
parameters P > P* allow to solve the problem (obviously, the parameter P* is 
problem-dependent). The problem 'How to determine P* ?' is not discussed in 
such an analysis, and in every concrete case it is solved using additional informa- 
tion about the problem. For example, in methods using a given Lipschitz constant 
(see survey |[T2l) P* is the Lipschitz constant and it is not discussed how to obtain 
it. Other examples are diagonal methods (see ||25| ) and information algorithms (see 
Il3n[38l ) using similar parameters. For all those methods it is possible to prove that 
if a value P > P* is used as the parameter, then they converge only to global min- 
imizers. An alternative approach is represented by methods converging to every 
point in the search domain (see |[T6ll44l '). 

In this paper a constrained Lipschitz global minimization problem is consid- 
ered. In an informal way it can be stated as follows: 

(i) The objective function is multiextremal, non-differentiable, 'black box', and 
requires a high time to be evaluated; 

(ii) Constraints are non-convex (or even multiextremal) and non-differentiable. 
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leading to a complex feasible region consisting of disjoint, non-convex sub- 
regions; 

(iii) Both the objective function and constraints are Lipschitz functions with un- 
known Lipschitz constants; 

(iv) Both the objective function and constraints may be partially defined, i.e., if a 
constraint is not satisfied at a point, the rest of constraints and the objective 
function may be not defined at that point. 

It can be seen from this statement that the problem belongs to the class of prob- 
lems considered in (37], thus statement (S2) will be used hereinafter. One promis- 
ing way to face problem (i)-(iv) is the information approach introduced by Stron- 
gin in |l38l[39lEl. It uses Peano type space-filling curves (see ll2l 151 l26l l38l l42]| 
for examples of usage of space-filling curves in mathematical programming) to re- 
duce the original Lipschitz multi-dimensional problem to a Holder univariate one 
(a comprehensive presentation of this approach can be found in |42|). Global opti- 
mization of Holder functions (see |l9l[l9l|45l) has given new tools for solving the 
reduced one-dimensional problem. Peano curves avoid constructions of support (or 
auxiliary) functions usually used in the multi-dimensional Lipschitz optimization 
(see, for example, |[T5][T7l[T8ll25]|27] |42 l and references given therein). Of course, 
if the user knows that the objective function is differentiable and/or the problem 
is convex, there is no sense to work with Peano curves and specific methods ex- 
plicitly using information about differentiability or convexity can be applied. In 
contrast, when you deal with non-differentiable 'black box' multiextremal prob- 
lems, it is not possible to work with sophisticated techniques using derivatives (or 
other strong a priori suppositions) and such a reduction to one dimension can help 
significantly. 

In this paper, a novel algorithm belonging to the family of information methods 
is proposed. It uses two powerful ideas for solving problem (i)-(iv). The first 
one is the index scheme (see ll35l [38l [39l 1411 l42l ). allowing to solve Lipschitz 
problems where both the objective function ip{y) and constraints Gi{y), I < i < 
m, may be multiextremal and partially defined. Its importance increases in this 
case because it is not clear how to solve such problems by using, for example, the 
penalty approach. In fact, the latter requires that ip{y) and Gi{y), 1 < i < m, are 
defined over the whole search domain. It seems that missing values can be simply 
filled in with either a big number or the function value at the nearest feasible point. 
Unfortunately, in the context of Lipschitz algorithms, incorporating such ideas can 
lead to infinitely high Lipschitz constants, causing degeneration of the methods and 
non-applicability of the penalty approach. Thus, for problems where the Lipschitz 
condition is almost a unique additional infonnation, the ability of the index scheme 
to work with partially defined problems becomes crucial. Moreover, the index 
scheme does not introduce any additional parameter and/or variable, whereas the 
penalty approach requiring determination of the penalty coefficient. 
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The second idea being at the basis of the new method is the local tuning on the 
behaviour of the objective function and constraints (see Il3n[32ll34l l35l). Original 
infoiTnation methods work with global adaptive estimates of the Lipschitz con- 
stants, i.e., the same estimates are used over the whole search region for ip{y) and 
functions 1 < i < m. However, global estimates (adaptive or given a 

priori) of the Lipschitz constants may provide a poor information about the be- 
havior of the objective function over every small sub-region of D. It has been 
shown in Il25ll3ni32l[34l[35l . for different classes of global optimization problems, 
that local estimates for different sub-regions of D can accelerate the search signif- 
icantly. Of course, it is necessary to balance global and local information about 
ip{y) obtained by the method during the search. Such a balancing is very important 
because using only the local information can lead to missing the global solution 
(seeEaiBll). 

The next section contains a formal statement of the problem (i)-(iv) and presents 
the new algorithm. Convergence conditions of the new method are established in 
Section 3. Numerical experiments collected in Section 4 show a satisfactory per- 
formance of the new algorithm in comparison with two global optimization tech- 
niques. Finally, Section 5 gives concluding remarks. 

2 Theoretical background and the index information 
algorithm with local tuning 

We start by formulating the Lipschitz optimization problem satisfying require- 
ments (i)-(iii). Find the constrained global minimum (p* and at least one minimizer 
y* such that 

if* = ifiy*) = minMy) : y G S, Gi{y) < 0, 1 < i < m}, (1) 
where the search domain is the hyperinterval 

S = {yGM^:a, <yj<6j, 1 < j < iV}, 

is the A^-dimensional Euclidean space, and the objective function ip{y) (hence- 
forth denoted as Gm+i(y)) and the functions Gi{y), 1 < i < m, of the con- 
straints can be multiextremal and non-differentiable, satisfying the Lipschitz con- 
dition with unknown constants < Lj < oo, 1 < i < m + 1, i.e., 

\G,{y')-Gi{y")\<U\\y' -y"\l l<i<m + l, y\y" e S. (2) 

Without loss of generality, we shall consider the search domain S = D, where 

= {y e : -2-1 < y, < 2-\ 1 < j < N}. (3) 

Formulation ([T]!-® assumes that all the functions Gi{y), I < i < m + 1, can be 
evaluated in the whole region D. In order to incorporate requirement (iv), we shall 



4 



assume that each function Gi (y) is defined and computable only in the correspond- 
ing subset Qi C D, where 

Qi = D, Qi+i = {y G : Gi{y) < 0}, l<i<m. (4) 

The above assumption also imposes the order in which the functions Gi{y), 1 < 
i < m, are evaluated. In many applications this order is determined by the nature 
of the problem. In other cases the user introduces a specific order suit for some 
reasons (for example, first verify easier computable constraints). In view of (01), 
the initial problem ([Til-© is rewritten as 

(p{y*) = mm{Gm+i{y) --y ^Qm+i}, (5) 

\Gi{y')-Gi{y")\<L,\\y' -y"\\, y',y"eQi, l<i<m + l, (6) 
where Li < Li. 

In order to start the description of the method it is necessary to remind the idea 
of the space-filling curves. Such curves were first introduced by Peano in f24] and 
Hilbert in [ 13] and are fractal objects constructed on the principle of self-similarity. 
They possess the property 'to fill' any cube D in M^, i.e., they pass through every 
point of D. An example of construction of such a curve (see ||29] l42l for details) 
is given in Fig. [T] Naturally, in numerical algorithms, approximations of the curve 
are used (Fig. [Upresents approximations of levels one to four). 

It has been shown in lH [38l |40l (see also 142 J ) that the multi-dimensional 
problem 

V^* = V{y*) = min{(/9(y) : y^D], (7) 

where Lp{y) is a Lipschitz function with constant L, < L < oo, can be reduced 
by a Peano-type space-filling curve y{x) to the one-dimensional problem 

r = fix*) = mini f{x), xe [0,1]}, (8) 

where the notation f{x) = ip{y{x)) is used for the obtained reduced one-dimensional 
function. Moreover, (see Theorem 8.1 in ||42] ). the function /(x) satisfies the 
Holder condition 

- f{x")\ <H\x'- x"\^"', x',x" e [0, 1], (9) 

in the Holder metric 

p{x',x") = \x' -x"\'/^, (10) 

where N is from © and = 2 L ^/N + 3. 

By applying the same curve y{x) to the functions Gi{y), 1 < i < m + 1, 
and using designations gi{x) = Gi{y{x)), the multi-dimensional problem Q-® 
is reduced to the following constrained one-dimensional problem (fTTI) -(fT3]): 

g*m+i = 9m+i{x*) = m.m.{gm+i{x) : X G Qm+l}, (11) 
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Figure 1 : Approximations of levels one to four to the Peano (Hilbert) curve in two 
dimensions. 



where the region Qm+i is defined by the relations 

qi = [0, 1], Qi+i = {x £ Qi-. gi{x) < 0}, l<i<m. (12) 

and the reduced functions gi{x), 1 < i < m + 1, satisfy the corresponding Holder 
conditions 

\giix')-gi{x")\<Hi\x' -x"\^/^, x',x"eqi, 1 < i < m + 1, (13) 

with the constants Hi = 2 Lj \/N + 3, where Lj,l<i<m + 1, are from (|6]l. 
Note that for the regions Qi from dUl and qi from ([T2l ) the following relation holds: 

Qi = {y{x) : X G gj, 1 < i < m + 1. 

This problem may be rewritten by using the index scheme proposed in |[39ll4T]| 
(see also |40|). The scheme is an alternative to traditional penalty methods. Instead 
of combining the objective and constraint functions into a penalty one (and the need 
to define in a proper way penalty coefficients), the scheme does not introduce any 
additional parameter or variable and evaluates constraints one at a time at every 
point where it has been decided to try to evaluate gm+i{x). The function gi{x) is 
calculated only if all inequalities 

gj{x)<0, l<j<i, 
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have been satisfied. In its turn, tiie objective function gm+i{x) is evaluated only 
for those points where all the constraints have been satisfied. 

The index scheme juxtaposes to every point of the interval [0, 1] an index 



z/ = z/(x), 1 < < m + 1, 



which is defined by the conditions 



gj{x) < 0, l<j<iy-l, gu{x) > 0, 



(14) 



where for = m + 1 the last inequality is omitted. Thus, u{x) represents the 
number of the first constraint not satisfied at x when u{x) < m + 1. Then = 
m + 1 means that all constraints were satisfied at x and the objective function 
gm+i{x) may be evaluated at this point. 

Let us now define an auxiliary function <I>(x) as follows 



where 5^+1 the solution of the problem (fTTI)-(fT3]). Due to (fT4l ) and (fTSl ). the 
function ^{x) has the following properties: 

i. <I>(x) > 0, when u{x) < m + 1; 

ii. <I>(x) > 0, when u{x) = m + 1; 

iii. ^{x) = 0, when z/(x) = m + 1 and gm+i{x) = g^+i' 

iv. <I>(x) is not continuous at the a priori unknown boundary points of the sets 
qi,2 < i < m + 1. 

Thus, the global minimizer of the constrained problem (ITTTi- ([T3l) coincides with 
the solution x* of the following unconstrained problem 



Of course, instead of the fractal Peano curve y{x) its d-level approximation yd{x) 
is used for evaluation of (approximations of levels one to four are shown in 
Fig-E from where it can be seen how a d-level approximation is obtained). 

The new index information global optimization algorithm with local tuning 
presented below generalizes and evolves two methods. On the one hand, the infor- 
mation global optimization algorithm with local tuning proposed in lISTl for solving 
the problem ^ and, consequently, the problem dH)-®. On the other hand, the in- 
dex algorithm with local tuning proposed in [35] for solving the one-dimensional 
problem 




z^(x) < m + 1 
z^(x) = m + 1 



(15) 



<^{x*) 



mm{^{x) : X € [0, 1]}. 



(16) 



min{/(x) : X £ [a,b], gi{x) <0, 1 < i < m} 
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where both f{x) and gi{x),l < i < m, are Lipschitz continuous functions. 

During every iteration of the algorithm a point x G [0, 1] is chosen and the 
value <I>(x) is evaluated (hereinafter such an evaluation will be called a trial and the 
corresponding point x a trial point). Suppose now that k iterations of the algorithm 
have already been executed (two initial trials are done at the end points x*^ = and 
x^ = 1). The point x^^^, A; > 1, is determined by the following algorithm. 

Step 1. The points . . . , of the previous iterations are renumbered by sub- 
scripts as follows 

= Xo < Xl < • • • < Xj < • • • < Xfc = 1. 

Step 2. To each point Xj associate the index Ui = v{xi), and the value 

Vi <m+l, 
i/j = m + 1, 

where 

zl = mhi{gm+i{xi) : < i < k,v{xi) = m + 1}. (17) 

The value estimates the unknown value g'^j^i from ([T5] ) on the basis of 
the available data. 

Step 3. Calculate lower bounds for the global Holder constants Hj of the func- 
tions (x) , 1 < J < m + 1 , as follows 

/ij = maxj— — ^iyr- : Q < q < p < k,Vp = Vq = j\ . (18) 
I- (xp - XgjVJv J 

Whenever /ij can not be calculated, set = 0. 

Step 4. For each interval [xj_i, Xj], i = 1, . . . , A;, calculate the values 

Mi = max{Ai,7i,e} (19) 

that estimate the local Holder constant over the interval [xi_i, Xj]. The val- 
ues Aj and 7^ reflect the influence on Mj of the local and global information 
obtained during the previous k iterations, ^ is a small number - parameter of 
the method. The values and 7j are defined by 

Ai = max{/i,Ci,ri}, (20) 
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where 



Ci = < 



0, otherwise 



(xi_i-a:,_2)V^ 

0, otherwise 



, if i > 2, Vi_2 = Vi-i-, > n 



if i < A; - 1, fj+i = i^i, z/j > 



and 



0, otherwise 



7j = > J = max{:/j, 



(21) 



where Xj"^ =m£ix{(a;i — Xj-i)-'^/^ : m£ix{z/j, fj-i} = j, 1 < ? < fe}. 

Step 5. For each interval Xj], i = 1, . . . , fc, calculate the characteristic of 

the interval 



Ai + 
2 A 



(zi-Zi-i)^ 2(zj+Zi_i) 



2A,; 



(rMi)2Ai 
_ ^ 

' rMi' 



rMi 



rMi 



(22) 



where Aj = (xj — .t,;_i)^/^ and r > 1 is a real value - the rehabihty 
parameter of the method. 

Step 6. Choose the interval [xt-i,xt] having the maximal characteristic as follows 



t = minjarg max{i?j . 1 <i < k}}. 

Step 7. If the interval [xt-i,xt] is such that 

{xt-xt-i)'/'' >S, 

where 5 is a given tolerance, go to Step 8. Otherwise Stop. 

Step 8. Execute the {k + l)-th iteration at the point 

fc+i ^ xt + xt-i _ sign{zt - zt-i)\zt - zt-i\^ 
^ ~ 2 2rMf 



(23) 



(24) 



if vt = ^t-i- In all other cases do it at the point 

x^+^={xt + xt-i)/2. (25) 

Set A; = A; + 1 and go to Step 1. 

If trial points with the index m + 1 have been generated by the algorithm then 
the value from ([TT] ) can be taken as an estimate of the global minimum 99* from 
([T]) and the corresponding point yd{x\) such that = Qm+ii^l.) as an estimate 
of the point y* . If no points with the index m + 1 have been generated by the 
algorithm then it is necessary to continue the search with changed parameters of 
the method. 

Let us give a few remarks on the algorithm introduced above. The information 
algorithms are derived as optimal statistical decision functions within the frame- 
work of a stochastic model representing the function to be optimized as a sample 
of a random function. The characteristic Ri in terms of the information approach 
(see ll38l l42l ) may be interpreted (after normalization) as the probability of find- 
ing a global minimizer within the interval based on the data available 
during the current iteration. The method uses in its work four parameters: d,r,^, 
and 6. Their choice wiU be discussed in Section 4 while presenting the numerical 
experiments. 

For every sub-interval Xj], I < i < k, global estimates fij of the global 

Holder constants Hj,l < j < m + 1, from ([T3T l are not used. In contrast, lo- 
cal estimates Mi,l < i < k, from ([T9] l are adaptively determined. The values 
Aj and 7^ reflect the influence on Mj of the local and global information obtained 
during the previous k iterations. When the interval is small, then 

is small too (see (|2TI )) and, due to ([T9l ), the local information represented by Aj 
has major importance. The value Aj is calculated by considering the intervals 
[xi-2, Xi-i], [xi-i,Xi], and [xj, Xj+i] (see (l20l )) as those which have the strongest 
influence on the local estimate. When the interval [xj_i, x^] is very wide, the local 
information is not reliable and the global information represented by 7j is used. 
Thus, local and global information are balanced in the values Mj ,1 < i < k. 
Note that the method uses the local information over the whole search region [0, 1] 
(and, consequently, over the whole multi-dimensional domain D) during the global 
search both for the objective function and constraints (being present in a implicit 
form in the auxiliary function ^>(x)). 

3 Convergence conditions 

In the further theoretical consideration it is assumed that each region Qi,2 < i < 
m + 1, from (|4ll is a union of Tj disjoint sub-regions Q^, 1 < i < Tj, having 
positive volume (in general, the numbers Tj are unknown to the user) where each 
sub-region Ql,l < j < Tj , 2 < f < m + 1 , is robust and can be non-convex. 
It is assumed that the feasible region Qm+i is nonempty. Recall that a set B is 
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robust if for each point j3 belonging to the boundary of B and for any e > there 
exists an e-neighborhood such that P| B has a positive volume. Suppo- 
sition about robustness of the sets Ql is quite natural in engineering applications 
(for example, in optimal control) where the optimal solution must have an admis- 
sible neighborhood of positive volume. This requirement is a consequence of the 
inevitable inaccuracy of the physical systems, and ensures that small changes of 
the parameters of the system will not lead the solution to leave the feasible region. 

Theoretical results presented in this section have the spirit of existence theo- 
rems. They show (similarly to the penalty approach) that for a given problem there 
exist parameters of the method allowing to solve this fixed problem. Recall that 
(see OTl ) the information available from the statement dS)-® is not sufficient for 
establishing concrete values of the parameters. These should be chosen from an ad- 
ditional information about the nature of the practical problem under consideration 
(see the next section). 

The first result links solutions of the original multi-dimensional problem 
(O to the solutions of the corresponding one-dimensional problem ([TTI) -([T3]) re- 
duced via the curve yd{x). It shows that for any given e > there exists a level d 
such that the approximation yd{x) passes through the e-neighborhood of the global 
solution y* of the problem ([S]!-®. Moreover, the global solution y*^ of the reduced 
problem ([TTI)-([T3]) will also belong to the same e-neighborhood. 

Theorem 1 For every problem (l?])-© having the unique solution y*, and a given 
accuracy e > 0, there exists a subdivision level d such that in the e-neighborhood 
e(y*) of the global minimizer y* there exists a segment h of the level-d approxi- 
mating curve yd{x) with the following properties: 

1. the segment h belongs to Qm+i H e(y*) and has a finite positive length; 

2. the segment h contains the global minimizer of the problem 

"PiVd) = va.m{Gm+i{yd{x)) : yd{x) G Qm+i,x G [0,1]}. (26) 

Proof. We have assumed that the feasible region Qm+i is the union of Tm+i 
robust subregions Q^+i, 1 < j < Tm+i- Suppose that y* £ Q^+v where Q^_^_l 
is one of those subregions. Then Q^^i H e{y*) is a robust set with a nonzero 
volume. This means that by increasing the accuracy of approximation it is possible 
to find an approximation level d such that the curve yd{x) will have a segment 
h C Qm+i ^ ^(y*) with the required properties. 

An illustration to Theorem [T]is given in Fig.|2] where a region Q^_|_i is shown 
in grey color Analogously, if problem dD-® has more than one solution, an 
approximation level d providing existence of such intervals having finite posi- 
tive lengths can be found for each of global minimizers. Of course, (l26l ) may 
be satisfied only for one interval but, anyway, all these intervals will contain e- 
approximations of the global solution. Thus, Theorem \T\ allows us to concentrate 
the further theoretical investigation on the behaviour of the method during the so- 
lution of problem (fT6l ). 
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Figure 2: Illustration of existence of a subdivision level d allowing an e- 
approximation of the global solution y* such that the point is the global 
minimizer of the objective function Gm+i{yd{x)) on Qm+i H yd{x). 

Lemma 1 Let x be a limit point of the infinite sequence {x^} generated by the 
algorithm and let i = i(k) be the number of an interval [xi^i,Xi\ containing this 
point during the k-th iteration. Then 

lim (xif^k) - Xi(k)-i)^^^ = 0, (27) 

fc— i-oo ^ ' ^ ' 

and for every o" > there exists an iteration number n{a), such that the inequality 

R^(k) < ^ (28) 

holds for all k > n{a). 

Proof. The new trial point x^^^ from (l24l) or (l25T l falls into an interval [xi-i, Xi] 
(where i{k) = t{k) is determined by (l23Tl). and divides this one into two sub- 
intervals 

Let us show that there exists a number a, independent of the iteration number k, 
such that 

max{xi - - Xj.i} < a{xi - Xj_i), 0.5 < a < 1, (29) 

holds for these intervals. In the case ut ^ ut-i the estimate (l29l ) is obtained im- 
mediately from (l25T l by taking ai = 0.5. In the opposite case, from ([T9l)-(l2TI) it 
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follows that 

<Xt- Xt-1, Vt = T^t-l- (30) 

From this inequality, (l24l) . and the fact that r > 1 we can conclude that in the case 
ut = ut-i ( |29l ) is true for a2 = (r + l)/(2r). Since 02 > oti, ( l29l ) is proved by 
taking a = 02- 

Now, considering (|29l ) together with the existence of a sequence converging to 
X (this point is a limit point of {x^}), we can deduce that (l27l) holds. Note that in 
the case when two intervals containing the point x there exist (i.e., when x G {x^}) 
the number i = i{k) is juxtaposed to the interval for which (|27] ) takes place. 

Let us prove (|28] ). From (l22l ) and (|30] | we obtain, for intervals with indexes 



Ri<{l+ r-^){xi - < 2{xi - Xi_i) 



l/Af 



taking into account that ^{x) > and r > 1. The last estimate holds also for 
intervals with Ui 7^ i^i-i- This inequality and (l27l) lead to 



Theorem 2 Lef x* fee any solution to the problem di6D awJ j = j (A;) fee the number 
of an interval [xj-i,Xj\ containing this point during the k-th iteration. Then, if for 
k > k* the condition 



holds, the point x* will be a limit point of {x'^}. The values Cj and Dj used in liSli 
are 



\ 0, otherwise. 
Proof. Consider the case Uj^i = uj. Due to (l32l) . we can write 

Zj.i<Cj{x* -Xj^i)'/^, (34) 

Zj <Cjixj -x*)^/^. (35) 
Now, by using (l34l) . (|35] ). and the designation 
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we deduce 

= + (1 _ - 

< Cj{xj - Xj^iY'^ max{/3V^ + (1 - p)^^ : < /? < 1} 
= 2i-V^C,(x,--x,_i)V^. 

By using this estimate and (|3TI ). (l33l) . we obtain 



{rMjYixj-Xj^iYl^ rMj 



> {xj-Xj.iY'^ {l+D'^ArMj)-^ -2'^-^/^ Cj{rMj)-^). 



Now, due to (|3TI ). we can conclude that 



> 0. (36) 



In the case Vj^i > Uj the estimate (1341 ) takes place due to (I321 ). From (I22l i it foUows 
that 

> 2{xj - x,_i)^/^(l - 2C,-(rM,-)-^), 

and, consequently, by taking into consideration (|3TI ). the inequality (l36l ) holds in 
this case also. Truth of (l36l ) for z^j_i < Vj is demonstrated by analogy. 

Assume now that x* is not a limit point of the sequence {x^}. Then, there 
exists a number Kg such that for all k > Kg the interval [xj-i,Xj\, j = j{k), is 
not changed, i.e., new points will not fall into this interval. 

Consider again the interval from Lemma [T] containing a limit point 



X. It follows from (1281 ) that there exists a number Kp such that 

for all k > k* = max{Ks, Kp}. This means that, starting from the iteration k*, 
the characteristic of the interval [xj_i, Xj], i = i{k), k > k*, is not maximal. Thus, 
a trial will fall into the interval [xj-i,Xj]. But this fact contradicts our assumption 
that X* is not a limit point. 

Corollary 1 Given the conditions of Theorem |2] Theorem [7] ensures that trial 
points generated by the algorithm will fall in the feasible sub-region Qm+i ne(y*). 

Theorems [T]|2] generalize for the constrained case (ITTt -(IT3Ti results established 
in f31] for the Lipschitz global optimization problems with box constraints. Usu- 
ally, Lipschitz global optimization algorithms need an overestimate (adaptive or a 
priori given) of the global Lipschitz constant for the whole search region (see, for 
example, survey |[T2l and references given therein). The new algorithm does not 
need knowledge of the precise Lipschitz constants for the functions Gi{y), 1 < 
i < m + I, over the whole search region. In fact, the point y* (as all points in D) 
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has up to 2 images on the Peano curve. To obtain an e-approximation of y* it is 
sufficient the fulfillment of the condition (|3TI ) (which is significantly weaker than 
the Lipschitz condition) for only one image of y* on a segment h of the curve such 
that/i C Qm+ine(y*). In the rest of the region D the constants Lj, 1 <i< m+1, 
can be underestimated. 

In the preceding it has been established that, if condition (|3T]) is satisfied, the 
trial sequence {x''} generated by the algorithm converges to the global minimizer 
X* of the function <&(y) and, consequently, the corresponding sequence {yd{x^)} 
converges to the global minimizer y* of the function 93(2/). In the following theo- 
rem it is proved that for any problem (IS)-® there exists a continuum of values of 
r satisfying (|3TI ). 

Theorem 3 For any problem (El-® ci value r* exists such that, for all r > r*, 
condition dJil ) is satisfied for all iteration numbers k >1. 

Proof. Let us choose 

r* = 23-i/A'^-i^jv + 3 max{Li : 1 < i < m + 1}, (37) 

and take a value r > r*. Since, due to ([T3l) . 

Hi = 2ViV + 3Lj, l<i<m + l, 

the value r* from (|37] | can be rewritten as 

r* =22"W^-i/7, 

H = maxjffj : 1 < i < m + 1}. 

For all k>l, for the estimates Mj(j.-^ and the values Cj^^^ from (|32] | the inequali- 
ties 

i < Mj{k)-, Cj(^k) < H, 
hold. Then, when the interval [xj-i,Xj] has Vj = uj^i^it follows from (|37] ) that 

To complete the proof it is sufficient to note that in the case Vj 7^ Uj^i the last 
estimate can be substituted by 

22-i/A^Cj >2Cj, 

because N > 2. Thus, for the chosen r > r* , condition (|3T]) is satisfied. 

Theorem |2] establishes sufficient convergence conditions for the method. The- 
orem |3]ensures that there exists a continuum of values of the parameter r satisfying 
these conditions. However, these theorems just prove existence of such values and 
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can not be considered as an instrument for determining the value r*. This value is 
problem-dependent and can not be found without additional information about the 
objective function and constraints. In presence of such information it is possible to 
solve the problem (SI), otherwise only the problem (S2) can be faced. 

Condition (|3T]) gives us a suggestion how to choose a reasonable value of r to 
start optimization in the case (S2). Note that the value Dj from (l34l ) can be equal 
to zero. The role of Mj is to estimate Cj. Thus, if Dj = and Mj = Cj, it 
follows from (OTI ) that r should be greater than 2^~^/^. Some practical advises for 
the choice of the method parameters will be given in the next section. 

4 Numerical experiments 

This section presents numerical experiments that investigate the performance of 
the proposed algorithm in solving some test problems. The first three problems are 
taken from the literature; the other problems are proposed by the authors. All the 
methods have been implemented in MATLAB 1 43 1 and the experiments have been 
executed at a PC with Pentium III 733MHz processor. 

Experiment 1. In the first experiment the new algorithm is compared with two 
methods: (i) the information algorithm for solving problems with box constraints 
proposed in |391, combined with the penalty approach; (ii) the original index in- 
formation algorithm proposed by Strongin and Markin in [39 1-|41 1, that does not 
use the local tuning. These methods have been chosen for comparison because all 
of them have similar computational cost for a single iteration (of course, if the cost 
required for the search of the penalty coefficient allowing to solve the problem is 
not taken in consideration for the method (i)), they use Peano curves, and have 
the same stopping rule. Since the penalty approach needs the objective function 
and constraints defined over the whole search region, test problems enjoying this 
property have been chosen. 

Problem 1. The first problem (see [41]) is to minimize the function 

ifiy) = -1.5 yl exp(l - - 20.25 (yi - ys)') 

- (0.5 (yi - l)(y2 - 1))^ exp(2 - (0.5 (yi - 1))^ - (ya - 1)^) 

(38) 

over the rectangle D = [0, 4] x [—1,3], under the constraints 

Gi(y) = 0.01 ((yi - 2.2)^ + (y^ - 1.2)2 _ 2.25) < 0, 
G2(y) = 100 (1 - (yi - 2)2/1.44 - (0.5 y2)2) < 0, 
G^iy) = 10 (y2 - 1.5 - 1.5 sin(6.283 (yi - 1.75))) < 0. 

The feasible region (see Fig. |5]l is the intersection of the areas inside the circle 
Gi{y) = 0, outside the ellipses G2(y) = 0, and below the sinusoide G3(y) = 0. 
It consists of three disjoint and non-convex pieces with non-smooth boundaries 
shown in Fig.|5]by the gray color. Here the solution is y* = (0.942, 0.944), giving 
the value ip{y*) = -1.489 |^. 



16 




Problem 2. In this problem (see BOll ) the same function as in the first experi- 
ment, over the same domain D, is minimized subject to the constraints 

Gi{y) = -{yi - 2.2)2 _ _ ^ < -1.21, 
G2(y) = (yi- 2.2)2 + (y2- 1.2)2 < 1.25, 

that define a narrow annulus centered at the point C = (2.2, 1.2). Here the solution 
is y* = (1.088, 1.088), giving the value ip{y*) = -1.477 Il42ll . 

Problem 3. This is a minimization problem with four constraints (see Section 
9.2 in [il4.n : the expression of the objective function is reported in the appendix; 
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0.5 1 1.5 2 2.5 3 3.5 4 

Figure 4: Trial points generated by the Strongin-Markin algorithm with 6 = 
for Problem 1 



the search region is D = [0, 80] x [0, 80] and the constraints are 

Gi(y) = 450 -yiy2<0, 

G2(y) = (0.1yi-l)2-y2<0, 

Gsiy) = 8 (yi - 40) - (ya - 30)(y2 - 55) < 0, 

Giiy) = {yi - 35) (yi - 30)/125 + 2/2 - 80 < 0. 

The feasible region for this problem is connected but not convex and has a non- 
smooth boundary. The solution is the point y* = (77.19, 64.06), and the value of 
the function at this point is ip{y*) = -59.59 [42|. 

Problem 4. The fourth test problem, proposed by the authors, is defined by the 
search domain D = [0, 2 vr] x [0, 2 tt], the non-smooth objective function 

ip{y) = -I sin(yi) sin(2 2/2)! + 0.01 (yi 2/2 + (2/1 - vr)^ + 3 (y2 - vr)^) , (41) 
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and the non-smooth constraints 

Giiy) = 1 - + vr/2 - I siii(2yi)| + yi/3 < 0, 

G2{y) = ?/2 - 37r/2 + 4 |(sin(yi + 7r)| + yi/3 - 1.9 < 0. 

The feasible region is made of two disjoint pieces, each of which is non-convex 
and has a non-smooth boundary. The best solution, found by a fine grid over the 
search region, is y* = (1.247, 2.392), yielding ip{y*) = -0.864. 



In all the experiments the maximum number of iterations was set to 5000, the 
level of approximation of the Peano curve to d = 10, the value ^ from Step 4 to 
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Table 1: Results for Experiment 1 with 6 = 10 ^ 







Strong 


lin-Markin method 








New method 






Problem 


711 


n2 


"-3 


714 






Til 


n2 




714 


71^ 




1 


1712 


953 


673 




279 


-1.489 


289 


201 


119 




60 


-1.455 


2 


1631 


672 






246 


-1.477 


316 


142 






31 


-1.438 


3 


351 


348 


316 


140 


88 


-59.59 


106 


103 


90 


65 


36 


-59.34 


4 


771 


359 






163 


-0.864 


115 


60 






29 


-0.863 



Table 2: Results for Experiments 1 with 6 = 10 







Strong 


;in-Markin alj 


;orithm 








New method 






Problem 


Til 


712 




714 


71^ 




711 


n2 


713 


714 


71,^ 




1 


4494 


2636 


1487 




826 


-1.489 


362 


274 


162 




102 


-1.489 


2 


4926 


1941 






781 


-1.478 


376 


170 






44 


-1.439 


3 


1073 


1070 


1032 


655 


359 


-59.59 


115 


112 


99 


74 


45 


-59.34 


4 


1867 


787 






355 


-0.864 


135 


72 






40 


-0.863 



^ = 10~ , and the value of the reliability parameter to r = 2.2. Only the points 
and 1 have been used as starting points. 

In the experiments with the penalty approach the constrained problems were 
reduced to the unconstrained ones as follows 

ifpiy) = ifiy) + P • max {Gi(y), G2{y), • • • , 0} ,y e D. 

The coefficient P has been computed by the rules: 

- the coefficient P has been chosen equal to 0.1 for all the problems and it has 

been checked whether the found solution for each problem belongs or not to 
the feasible subregions; 

- if it does not belong to the feasible subregions, the coefficient P* has been 

iteratively increased by 0.1 until a feasible solution has been found. 

Figures [3]-l5] present trial points generated by the three methods during Exper- 
iment 1 with 6 = 10^^. Figure [3] shows the trial points generated by the penalty 
approach. Remind that both the objective function and all the constraints have been 
evaluated at every point. Figures |4] and |5] show the points where the algorithms us- 
ing the index approach have chosen to evaluate the constraints and eventually the 
target function. Points where the target function has been computed are denoted by 
stars; those where only the first constraint has been evaluated (and found violated) 
are denoted by crosses; circles denote points where the first constraint was satisfied 
and the second was not; squares denote points were the first two constraints were 
satisfied and the third was not. 
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Table 3: Results obtained by the information algorithm combined with the penalty 
approach Experiment 1 



Test Problem 


P* 


S = 


10"'' 


S = 


10-* 






n 




n 




1 


0.1 


2298 


-1.489 


5000 


-1.489 


2 


0.3 


428 


-1.478 


1621 


-1.478 


3 


0.8 


5000 


-59.59 


5000 


-59.59 


4 


0.5 


283 


-0.864 


789 


-0.864 



Table 4: Speed-up obtained in Experiment 1 for Problems 1-4 



Problem 


1 


2 


3 


4 


1 


2 


3 


4 






5 = 


IQ-^ 






5 = 


10-^ 




Si 


4.29 


5.13 


3.74 


4.03 


12.25 


10.84 


14.68 


8.83 




4.65 


7.94 


2.44 


5.62 


8.10 


17.75 


7.98 


8.88 


S3 


5.41 


5.21 


3.11 


6.34 


10.49 


12.96 


9.41 


11.27 


Si 


7.95 


1.35 


47.17 


2.46 


13.81 


4.31 


43.48 


5.84 


§2 


38.30 


13.81 


138.89 


9.76 


49.02 


36.84 


111.11 


19.73 


S3 


13.74 


2.63 


62.50 


4.16 


22.47 


8.24 


56.18 


9.58 



Table [T] compares performance of the proposed algorithm to the Strongin- 
Markin method for S = 10^^. The data reported are the computed solution ipc, 
the number of evaluations of each constraint rij, and of the target function n^. To 
investigate the sensitivity of both algorithms to the value of 5 the tests have been 
repeated with S = 10"^; results are reported in Table |2] Finally, Table |3]presents 
results obtained by the information algorithm combined with the penalty approach 
on Problems 1 — 4, where n is the number of iterations executed by the algorithm. 
Cases where the algorithm was not able to stop in 5000 iterations are shown in 
bold. 

Table 0] illustrates acceleration reached by the new algorithm in comparison 
with the two methods tested in Experiment 1 . Three speed-up indexes have been 
used: the first index is the ratio Si of iterations ni of the Strongin-Markin method 
over the new one; the second is the ratio 52 of of the objective function eval- 
uations; finally, S3 is the ratio of the sums n^p + YlT=i ^fc' being the number of 
summary evaluations of the objective function and constraints. The same quanti- 
ties for the penalty approach are indicated hy Si,l < i < 3. It can be seen from 
Table |4]that the new method significantly outperforms both methods tested. More- 
over, the computational burden increases very slowly for decreasing 6 for it, while 
the Strongin-Markin algorithm and the penalty approach show a fast increase of 
iterations when a better accuracy is required. 

Since Table|4]shows an evident superiority of the new method over the penalty 
algorithm, the latter is not used in the further experiments (also because it is diffi- 
cult to determine the correct penalty coefficients for problems used in these exper- 
iments). 
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Table 5: Results for Experiment 2 





r = 2.5 




r = 3 






Strongin-Markin method 


New method 


Strongin-Markin method 


New method 


Successes 


19/20 


18/20 


20/20 


20/20 


ni average 


2177 


359 


2939 


651 


n2 average 


874 


179 


1150 


330 


riip average 


282 


46 


312 


88 


Si 




6.06 




4.51 


S2 




6.13 




3.55 


S3 




5.71 




4.12 



Experiment 2. This experiment aims at investigating the sensitivity of both 
algorithms to the value of the reliability parameter r. We have generated twenty 
pairs of random numbers in the range (—1, 1), and have added any of them to the 
center C = (2.2, 1.2) of the annulus of Problem 2, thus obtaining twenty new 
problems, whose solutions have been found with a fine grid over the domain D. 

They have been solved with both methods for either r = 2.5 and r = 3. With 
r = 2.5 the proposed algorithm has found the solution in 18 cases out of 20, the 
Strongin-Markin algorithm was successful in 19 cases out of 20. Note that the pro- 
posed algorithm has found the solution in the case the Strongin-Markin algorithm 
has failed. With r = 3, both algorithms have found the solution in all cases. The 
results (number of successes, average number of computation of the constraints 
and of the objective function) are reported in Table |5] With r = 3, in one case 
the Strongin-Markin algorithm has reached the maximum number of iterations, but 
anyway it has found a good approximation to the solution. It can be seen from Ta- 
ble [5] that the computational burden increases very slowly for increasing r for the 
proposed algorithm, while the Strongin-Markin algorithm shows a fast increase of 
iterations. 



Experiment 3. The last experiment involves a problem having changeable 
dimension N and three constraints. The objective function is 

N 

VN{y) =yi + eP--\''"^-'"'-+^\, pN = (Y^ yf)'^. (42) 

j=i 

The search region is D = [—2, 8] x [—6, 4] x . . . x [—6, 4] and the constraints are 

Gi (y) = (2/1 - 2.5)2 - 6.25 + 2/' < 0, 
G2(y) = -(yi- 2)2 + 2.25 -EL 2/2 

Gsiy) = 2/2 - 37r/2 + 4 |(sin(yi + TT)\+yi/3 + 0.8 < 0. 

Numerical experiments have been executed with 2 < < 6 with different values 
of r and 6. The level of approximation of the Peano curve was set to d = 10 for 
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2 < < 5 and to (i = 8 for = 6, the value ^ = 10 . In all the experiments 
the maximum number of iterations was set to 40000. The global solution is y* = 
(0, . . . , 0), yielding ifNiv*) = exp(-4) ~ 0.0183. 

The obtained results are shown in Tables |6] and |7] In all the experiments the 
new algorithm is significantly faster. It can be seen from the Table |7] that starting 
from dimension = 3 the Strongin-Markin algorithm was not able to stop in less 
than 40000 iterations. 

Let us comment the choice of parameters of the new method. As it has been 
mentioned in the Introduction, the statement of the problem does not allow to create 
any deterministic or stochastic algorithm that, after a finite number of evaluations 
of ip{y), is able to provide an underestimate of the global minimum (p{y*). These 
means that it is impossible to say without some additional information about the 
problem which values of parameters ensure that the method finds the global mini- 
mizes However, the executed experiments can advice some recommendations. 

First of all, the parameters d and ^ are quite technical and can be chosen easily. 
To have a better accuracy of the solution it is necessary to have a better approxima- 
tion of the Peano curve determined by d. Since two points in the multidimensional 
space have two images at the interval [0, 1], then in a computer realization (see, 
for example, [42]) for a correct work of the method these images should be repre- 
sented by two different numbers. The mapping is such that for a dimension A^ and 
a level d two numbers a,b £ [0, 1] will be considered different if |a — 6| > 2"'^^ . 
For example, in double precision the minimal representable number is 2^^^. Thus, 
if one works with a simple realization of the Peano curve using double precision, 
the product dN should be less than 52. Of course, more sophisticated implementa- 
tions of the Peano curve realizing numbers with more digits allow to increase this 
number. 

The second parameter is ^ and it is chosen equal to a small number (see experi- 
ments). It has been introduced in the method to ensure that it works correctly in the 
case li = Ci = ri = where k, Ci,ri are from (|20] |. The case when the real Holder 
constants are less than ^ is degenerous for the method because the local tuning is 
not used and the method works slower. 

The parameters r and 6 are chosen by the user on the basis of additional in- 
formation about the problem or just by increasing r and decreasing 6 (see Table [6l) 
trying to obtain a satisfactory result (see statement (S2) in Introduction). It can be 
also seen from Table [6]that by increasing dimension of the problem it is necessary 
to change r and 6 in the same way. This happens because, in spite of increasing 
dimension, the reduced problem is always determined over the interval [0, 1] inde- 
pendently on the dimension A^. As a consequence, the reduced problem has more 
local minima, and to find the global one it is necessary to increase reliability of the 
method (by increasing r) and to make the search more accurate (by decreasing 6). 
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Table 6: Results for the new method on Problems (|42] | with dimension two to six 



r 




S 


yi 


y2 


y3 


y4 




ya 


m 


712 


ns 






2.35 




-'S 


0.0034 


0.0596 










130 


51 


41 


36 


0.0295 


2.35 


10" 


-4 


0.0007 


0.0596 










147 


57 


47 


42 


0.0268 


2.45 


10" 


-3 


0.0068 


0.1185 


0.1182 








819 


470 


377 


249 


0.0555 


2.45 


10" 


-4 


0.0052 


-0.0088 


0.0107 








5847 


4723 


3975 


3313 


0.0252 


2.7 


10" 


~3 


0.2540 


-0.1553 


-0.0674 


0.4600 






708 


533 


364 


249 


0.6247 


2.7 


10" 


-4 


0.0264 


-0.0771 


0.0400 


-0.0679 






14441 


12485 


10081 


7410 


0.0621 


3.3 


10" 


-3 


0.5146 


-0.1650 


-0.0771 


-0.0924 


-0.3799 




620 


206 


164 


111 


1.1704 


3.3 


10" 


-4 


0.0654 


0.0303 


-0.2920 


0.0596 


-0.0692 




18535 


12289 


10761 


7344 


0.1748 


3.35 


10" 


-3 


0.5977 


-0.1992 


-0.0039 


-0.3298 


1.0898 


-0.6289 


445 


39 


35 


25 


1.9579 


3.35 


10" 


-4 


0.3193 


0.1523 


-0.3945 


0.1523 


0.3477 


0.1523 


6610 


2486 


2292 


1534 


0.9688 


3.35 


5 • 10" 


-5 


0.0508 


-0.0430 


0.0742 


-0.0430 


-0.1992 


0.0504 


19105 


9871 


8989 


6233 


0.1208 



Table 7: Results for the Strongin-Markin algorithm on Problems (l42l ) with dimen- 



sion two to six 



r 




S 


yi 


y2 


2/3 


y4 


2/5 


ya 


m 


712 


713 






2.35 


10" 


-S 


0.0025 


0.0010 










5071 


3811 


3314 


3066 


0.0211 


2.35 


10" 


-4 


0.0000 


0.0010 










39883 


33277 


31108 


30423 


0.0185 


2.45 


10" 


-3 


0.0388 


0.0010 


0.0498 








11003 


5734 


4911 


3736 


0.0654 


2.45 


10" 


-4 


0.0068 


0.0400 


-0.0351 








40000 


22591 


19901 


16023 


0.0320 


2.7 


10" 


-3 


0.1729 


0.0693 


-0.1162 


-0.1860 






5067 


1833 


1680 


1037 


0.2676 


2.7 


10" 


~4 


0.0459 


0.1210 


-0.2139 


0.0693 






40000 


17699 


15402 


10566 


0.1271 


3.3 


10" 


-3 


0.5635 


0.3037 


-0.0381 


0.3916 


0.0034 




2661 


432 


381 


220 


1.5087 


3.3 


10" 


-4 


0.1436 


0.3330 


-0.2627 


-0.0869 


0.0492 




40000 


10487 


9661 


6137 


0.3763 


3.35 


10" 


-3 


0.4756 


-0.2920 


-0.1846 


-0.0186 


-0.4482 


-1.1709 


650 


23 


22 


18 


1.9356 


3.35 


10" 


~4 


0.3877 


0.1377 


0.4111 


0.2061 


-0.3115 


-0.1260 


40000 


5413 


5091 


3251 


1.1450 


3.35 


5 ■ 10" 


~5 


0.3877 


0.1377 


0.4111 


0.2061 


-0.3115 


-0.1260 


40000 


5413 


5091 


3251 


1.1450 
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5 Conclusions 



In this paper, a novel global optimization algorithm for solving multi-dimensional 
Lipschitz global optimization problems with multiextremal partially defined con- 
straints and objective function has been presented. The new algorithm uses Peano 
type space-filling curves and the index scheme to reduce the original constrained 
problem to a Holder one-dimensional problem. Local tuning on the behaviour of 
the objective function and constraints is executed during the work of the global 
optimization procedure in order to accelerate the search. The new method works 
without introducing any penalty coefficients and/or auxiliary variables. Conver- 
gence conditions of a new type have been established for the algorithm. 
The new algorithm enjoys the following properties: 

- in order to guarantee convergence to a global minimizer y* it is not nec- 
essary to know the precise Lipschitz constants for the objective function and 
constraints for the whole search region; on the contrary, only fulfillment of 
condition (|3TI ) for a segment of the Peano curve containing one of the images 
of y* is required; 

- the usual problem of determining the moment to stop the global procedure 
in order to start a local search does not arise because local information is 
taken into consideration throughout the duration of the global search; 

- local information is taken into account not only in the neighborhood of 
a global minimizer but also in the whole search region, thus permitting a 
significant acceleration of the search; 

- thanks to usage of the index scheme, the new algorithm does not introduce 
any additional parameters and/or variables. Constraints are evaluated at ev- 
ery point one at a time until the first violation of one of them, after that the 
rest of constraints and the objective function are not evaluated at this point. 
In its turn, the objective function is evaluated only for that points where all 
the constraints have been satisfied. 

The algorithm has been tested on a number of problems taken from literature. 
Numerical results show the good performance of the new technique in comparison 
with the two methods taken from literature. 

Appendix 

The target function for Problem 3 is 

^{y) = -{Bi + B2 yi + S3 yf + B4yf + B5 yf + ^6^2 + Bj yiy2 + Bg yly2 
+ Bg yfy2 + Bio yty2 + Bu yl + Bu yl + B13 y| + Bii/{1 + 7/2) 
+ Bi5 ylyl + B16 ylyl + Bn yfy^ + Bis Vivl + -^19 yiyf 
+ S2oexp(0.0005yi?/2)) 
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where coefficients Bi,. . . , B20 are 

Bi = 75.1963666677; B2 = -3.8112755343; B3 = 0.1269366345; B4 = -0.0020567665; 
B5 = 0.0000103450; Be = -6.8306567613; Br = 0.0302344793; Bs = -0.0012813448; 
Bq = 0.0000352559; Bw = -0.0000002266; Bn = 0.2564581253; B12 = -0.0034604030; 
Bi3 = 0.0000135139; Bu = -28.1064434908; B15 = -0.0000052375; Bu = -0.0000000063; 
S17O = 0.0000000007; Big = 0.0003405462; S19 = -0.0000016638; B20 = -2.8673112392. 
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